program dipole
implicit none

real :: x,y,z,V,V1,V2,dVx,dVy,Ex,Ey,E,ux,uy
integer ::i,j


do x=-20,20
  do y=-20,20
  V= potential(x,y)
  V1 = potential(x+.01,y)
  V2= potential(x,y+.01)
  dVx = V1-V
  dVy = V2-V
  Ex= -dVx/.01
  Ey= -dVy/.01
  E = sqrt(Ex**2+Ey**2)
  ux = Ex/E
  uy = Ey/E
  open(4,file='E_Plot')
  write(4,*) x,y,ux,uy

 end do
end do

contains
real function potential (x,y)
implicit none

real,intent(in)::x,y
real::px,py
px = 1
py = 0
potential = (px*x+py*y)/(sqrt(x**2+y**2)**3)

end function

end program

! command  for plot >>  plot 'filename' with vectors